Nonlinear Regression with R

setup

# devtools::install_github("cran/drc", force = TRUE)
# devtools::install_github("cran/alr3", force = TRUE)
# devtools::install_github("cran/NRAIA", force = TRUE)
# devtools::install_github("cran/nlrwr", force = TRUE)
library(nlrwr)

Loading required package: alr3
Loading required package: car
Loading required package: carData
Loading required package: drc
Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:alr3’:

    forbes


'drc' has been loaded.

Please cite R and 'drc' if used for a publication,
for references type 'citation()' and 'citation('drc')'.


Attaching package: ‘drc’

The following objects are masked from ‘package:stats’:

    gaussian, getInitial

Loading required package: HydroMe
Loading required package: lattice
Loading required package: lmtest
Loading required package: zoo

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

Loading required package: NISTnls
Loading required package: nlme
Loading required package: nls2
Loading required package: proto
Loading required package: nlstools

'nlstools' has been loaded.

IMPORTANT NOTICE: Most nonlinear regression models and data set examples
related to predictive microbiolgy have been moved to the package 'nlsMicrobio'

Loading required package: NRAIA
Registered S3 method overwritten by 'NRAIA':
  method           from 
  plot.profile.nls stats

Attaching package: ‘NRAIA’

The following object is masked from ‘package:nlstools’:

    plotfit

Loading required package: sandwich

'nlrwr' has been loaded.

Please cite R and 'nlrwr' if used for a publication,
for references type 'citation()' and 'citation('nlrwr')'.

Attaching package: ‘nlrwr’

The following object is masked from ‘package:nlstools’:

    confint2

1 Introduction

1.1 A stock-recruitment model

## Fig. 1.1. 
plot(num.fish ~ spawn.biomass, 
     data = M.merluccius, 
     xlab = "Spawning biomass (1000 tonnes)",
     ylab = "Recruitment (million fish)")

1.2 Competition between plant biotypes

## Fig. 1.2.
plot(biomass ~ x, 
     data = RScompetition, 
     log = "", 
     xlab = Density ~ (plants/m^2), 
     ylab = Biomass ~ of ~ sensitive ~ biotype ~ (g/plant), 
     pch = as.numeric(as.factor(RScompetition$z)))

1.3 Grouped dose-response data

## Fig. 1.3.
xyplot(DryMatter ~ Dose | Herbicide,
       data = S.alba, 
       scales = list(x = list(log = TRUE)),
       ylab = "Dry matter (g/pot)",
       xlab = "Dose (g/ha)")

2 Getting Started

2.2 Getting started with nls()

2.2.1 Introducing the data example

L.minor

## Fig. 2.1.
plot(rate ~ conc, data = L.minor,
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3))

2.2.2 Model fitting

L.minor.m1 <- nls(rate ~ Vm * conc/(K + conc), 
                  data = L.minor, 
                  start = list(K = 20, Vm = 120), 
                  trace = TRUE)

624.3282    (1.32e+00): par = (20 120)
244.5460    (2.24e-01): par = (15.92382 124.5715)
234.5198    (2.91e-02): par = (17.25299 126.4388)
234.3595    (5.72e-03): par = (17.04442 125.9618)
234.3533    (1.11e-03): par = (17.08574 126.0467)
234.3531    (2.15e-04): par = (17.07774 126.0302)
234.3531    (4.19e-05): par = (17.0793 126.0334)
234.3531    (8.14e-06): par = (17.07899 126.0328)

deviance(L.minor.m1)

[1] 234.3531

d2 <- sum((L.minor$rate -predict(L.minor.m1))^2)
all.equal(deviance(L.minor.m1), d2)

[1] TRUE

logLik(L.minor.m1)

'log Lik.' -24.86106 (df=3)

coef(L.minor.m1)

        K        Vm 
 17.07899 126.03276 

summary(L.minor.m1)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 8.141e-06

2.2.3 Prediction

fitted(L.minor.m1)

[1]  18.06066  28.56474  38.52679  71.09461  78.03806  87.78424  91.62683 116.28685
attr(,"label")
[1] "Fitted values"

concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 10))
predict(L.minor.m1, newdata = data.frame(conc = concVal))

 [1]  18.06066  75.09908  92.70509 101.26607 106.32776 109.67157 112.04518
 [8] 113.81733 115.19094 116.28685

2.2.4 Making plots

## Fig. 2.2.
plot(rate ~ conc, 
     data = L.minor,
     ylim = c(10, 130), 
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3)) 
lines(L.minor$conc, fitted(L.minor.m1))

## Fig. 2.3.
plot(rate ~ conc, 
     data = L.minor,
     ylim = c(10, 130), 
     ylab = "Uptake rate (weight/h)", 
     xlab = Substrate ~ concentration ~ (mmol ~ m^-3))

concVal <- with(L.minor, seq(min(conc), max(conc), length.out = 100))
lines(concVal, predict(L.minor.m1, newdata = data.frame(conc = concVal)))

abline(h = coef(L.minor.m1)[2], lty = 2)

2.2.5 Illustrating the estimation

L.minor.m1con <- nlsContourRSS(L.minor.m1)

0%100%
 RSS contour surface array returned 

## Fig. 2.4. 
par(pty = "s")
plot(L.minor.m1con, col = FALSE, nlev = 10)

2.3 Generalised linear models

## this code should be runned in the clean R session
## library(nlrwr) cause error 
# data(L.minor, package = "nlrwr")
# L.minor.m4 <- glm(rate ~ I(1/conc), 
#                   data = L.minor, 
#                   family = gaussian("inverse"))
# summary(L.minor.m4)

Excercises

## 2.1
y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)


e2.1 <- nls(RGR ~ y0 * exp(Day/b),
            data = RGRcurve,
            start = list(y0 = y0_start, b = b_start))

plot(RGRcurve$Day, RGRcurve$RGR, xlab = "Day", ylab = "RGR")
nd <- with(RGRcurve, seq(min(Day), max(Day), length.out = 101))
lines(nd, predict(e2.1, newdata = data.frame(Day = nd)), col = 2)

## 2.3
plotfit(L.minor.m1)

3 Starting Values and Self-starters

3.1 Finding starting values

3.1.1 Graphical exploration

## Fig. 3.1.
# install.packages("NISTnls")
data(Chwirut2, package = "NISTnls")
plot(y ~ x, 
     data = Chwirut2, 
     xlab = "Metal distance", 
     ylab = "Ultrasonic response")

## Fig. 3.2.
expFct <- function(x, beta1, beta2, beta3) {
  exp(-beta1 * x)/(beta2 + beta3 * x)
}

plot(y ~ x, 
     data = Chwirut2, 
     xlab = "Metal distance", 
     ylab = "Ultrasonic response",
     ylim = c(0, 100)
     )
curve(expFct(x, beta1 = 1, beta2 = 0.01, beta3 = 1), add = TRUE, col = 2)

## Fig. 3.3.
plot(y ~ x, data = Chwirut2, xlab = "Metal distance",
     ylab = "Ultrasonic response", ylim = c(0, 100))
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 1), add = TRUE, lty = 2,
      col = 2)

curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.1), add = TRUE, lty = 3,
      col = 3)
curve(expFct(x, beta1 = 0.1, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 4,
      col = 4)

curve(expFct(x, beta1 = 0.2, beta2 = 0.01, beta3 = 0.01), add = TRUE, lty = 1,
      col = 6)

Chwirut2.m1 <- nls(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = list(beta1 = 0.2, beta2 = 0.01, beta3 = 0.01)
)
summary(Chwirut2.m1)

Formula: y ~ expFct(x, beta1, beta2, beta3)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 0.1665753  0.0383032   4.349 6.56e-05 ***
beta2 0.0051653  0.0006662   7.753 3.54e-10 ***
beta3 0.0121501  0.0015304   7.939 1.81e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.172 on 51 degrees of freedom

Number of iterations to convergence: 5 
Achieved convergence tolerance: 4.684e-06

3.1.2 Searching a grid

## bad starting example
# Chwirut2.m0 <- nls(
#   y ~ expFct(x, beta1, beta2, beta3),
#   data = Chwirut2,
#   start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
# summary(Chwirut2.m0)

## nlsLM may solve the bad starting
Chwirut2.m0 <- minpack.lm::nlsLM(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = list(beta1 = 1, beta2 = 0.01, beta3 = 1))
summary(Chwirut2.m0)

Formula: y ~ expFct(x, beta1, beta2, beta3)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 0.1665765  0.0383033   4.349 6.56e-05 ***
beta2 0.0051653  0.0006662   7.753 3.54e-10 ***
beta3 0.0121500  0.0015304   7.939 1.81e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.172 on 51 degrees of freedom

Number of iterations to convergence: 22 
Achieved convergence tolerance: 1.49e-08

grid.Chwirut2 <- expand.grid(beta1 = seq(0.1, 1, by = 0.1),
                             beta2 =  0.01,
                             beta3 = seq(0.1, 1, by = 0.1))
grid.Chwirut2

Chwirut2.m2a <- nls2(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = grid.Chwirut2,
  algorithm = "brute-force"
)
Chwirut2.m2a

Nonlinear regression model
  model: y ~ expFct(x, beta1, beta2, beta3)
   data: Chwirut2
beta1 beta2 beta3 
 0.10  0.01  0.10 
 residual sum-of-squares: 60696

Number of iterations to convergence: 100 
Achieved convergence tolerance: NA

nls(
  y ~ expFct(x, beta1, beta2, beta3),
  data = Chwirut2,
  start = as.list(coef(Chwirut2.m2a)))  # starting value from Chwirut2.m2a

Nonlinear regression model
  model: y ~ expFct(x, beta1, beta2, beta3)
   data: Chwirut2
   beta1    beta2    beta3 
0.166578 0.005165 0.012150 
 residual sum-of-squares: 513

Number of iterations to convergence: 6 
Achieved convergence tolerance: 4.429e-06

3.2 Using self-starter functions

3.2.1 Built-in self-starter functions for nls()

L.minor.m2 <- nls(rate ~ SSmicmen(conc, Vm, K), data = L.minor)
summary(L.minor.m2)

Formula: rate ~ SSmicmen(conc, Vm, K)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
Vm  126.033      7.173  17.570 2.18e-06 ***
K    17.079      2.953   5.784  0.00117 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 2.478e-06

3.2.2 Defining a self-starter function for nls()

expModel <- function(predictor, b, y0) {
  y0 * exp(predictor/b)
}

expModelInit <- function(mCall, LHS, data) {
  xy <- sortedXyData(mCall[["predictor"]], LHS, data)
  lmFit <- lm(log(xy[, "y"]) ~ xy[, "x"])
  coefs <- coef(lmFit)
  y0 <- exp(coefs[1])
  b <- 1/coefs[2]
  value <- c(b, y0)
  names(value) <- mCall[c("b", "y0")]
  value
}

SSexp2 <- selfStart(expModel, expModelInit, c("b", "y0"))

with(RGRcurve, SSexp2(Day, 4, 0.2))

 [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2568051
 [8] 0.2568051 0.2568051 0.2568051 0.2568051 0.3297443 0.3297443 0.3297443
[15] 0.3297443 0.3297443 0.3297443 0.4234000 0.4234000 0.4234000 0.4234000
[22] 0.4234000 0.4234000 0.5436564 0.5436564 0.5436564 0.5436564 0.5436564
[29] 0.5436564 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378 0.8963378
[36] 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112 1.4778112

Warning: selfStart caused an ERROR

## ERROR
# getInitial(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)

#RGRcurve.m1 <- nls(RGR ~ SSexp2(Day, b, y0), data = RGRcurve)
#coef(RGRcurve.m1)

More examples on defining self-starter functions for nls() are found in Watkins and Venables (2006) or Venables and Ripley (2002a, pp. 216–217).

4 More on nls()

4.2 Supplying gradient information

4.2.1 Manual supply

MMfct1 <- function(conc, K, Vm) {
  numer <- Vm * conc
  denom <- K + conc
  mean <- numer/denom
  partialK <- -numer/(denom^2)
  partialVm <- mean/Vm
  attr(mean, "gradient") <- cbind(partialK, partialVm)
  return(mean)
}
MMfct1(1, 2, 3)

[1] 1
attr(,"gradient")
       partialK partialVm
[1,] -0.3333333 0.3333333

L.minor.mgr1 <- nls(rate ~ MMfct1(conc, K, Vm), 
                    data = L.minor, 
                    start = list(K = 1, Vm = 1),
                    trace = TRUE)

42634.85    (7.03e+00): par = (1 1)
31698.57    (9.47e+00): par = (389.6406 100.9883)
21357.28    (1.23e+01): par = (31.6085 46.75042)
18634.08    (1.03e+01): par = (27.26998 51.07311)
13989.24    (8.02e+00): par = (22.1864 59.6233)
7637.536    (5.64e+00): par = (18.35113 75.51722)
2003.136    (2.75e+00): par = (17.01927 100.4895)
234.3615    (6.31e-03): par = (17.1086 126.0567)
234.3533    (9.56e-04): par = (17.07326 126.021)
234.3531    (1.86e-04): par = (17.08017 126.0352)
234.3531    (3.61e-05): par = (17.07882 126.0324)
234.3531    (7.02e-06): par = (17.07909 126.0329)

summary(L.minor.mgr1)

Formula: rate ~ MMfct1(conc, K, Vm)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 11 
Achieved convergence tolerance: 7.025e-06

## without derivative
nls(rate ~ Vm * conc/(K + conc), 
    data = L.minor, start = list(K = 1, Vm = 1), 
    trace = TRUE)

42634.85    (7.03e+00): par = (1 1)
31698.57    (9.47e+00): par = (389.6406 100.9883)
21357.27    (1.23e+01): par = (31.60848 46.75041)
18634.07    (1.03e+01): par = (27.26996 51.07311)
13989.24    (8.02e+00): par = (22.18639 59.6233)
7637.536    (5.64e+00): par = (18.35113 75.51722)
2003.136    (2.75e+00): par = (17.01927 100.4895)
234.3615    (6.31e-03): par = (17.1086 126.0567)
234.3533    (9.56e-04): par = (17.07326 126.021)
234.3531    (1.86e-04): par = (17.08017 126.0352)
234.3531    (3.61e-05): par = (17.07882 126.0324)
234.3531    (7.03e-06): par = (17.07909 126.0329)
Nonlinear regression model
  model: rate ~ Vm * conc/(K + conc)
   data: L.minor
     K     Vm 
 17.08 126.03 
 residual sum-of-squares: 234.4

Number of iterations to convergence: 11 
Achieved convergence tolerance: 7.029e-06

4.2.2 Automatic supply

MMfct2 <- deriv(~ Vm * conc / (K + conc), 
                c("K", "Vm"), 
                function(conc, K, Vm) {})
MMfct2

function (conc, K, Vm) 
{
    .expr1 <- Vm * conc
    .expr2 <- K + conc
    .value <- .expr1/.expr2
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("K", 
        "Vm")))
    .grad[, "K"] <- -(.expr1/.expr2^2)
    .grad[, "Vm"] <- conc/.expr2
    attr(.value, "gradient") <- .grad
    .value
}

MMfct2(1, 2, 3)

[1] 1
attr(,"gradient")
              K        Vm
[1,] -0.3333333 0.3333333

L.minor.mgr2 <- nls(
  rate ~ MMfct2(conc, + K, Vm), 
  data = L.minor, 
  start = list(K = 20, Vm = 120))
summary(L.minor.mgr2)

Formula: rate ~ MMfct2(conc, +K, Vm)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
K    17.079      2.953   5.784  0.00117 ** 
Vm  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 8.148e-06

Ref: Venables와 Ripley(2002a, pp. 214–216)

4.3 Conditionally linear parameters

4.3.1 nls() using the “plinear” algorithm

L.minor.m3 <- nls(rate ~ conc/(K + conc), 
                  data = L.minor, 
                  algorithm = "plinear",
                  start = list(K = 20))
summary(L.minor.m3)

Formula: rate ~ conc/(K + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K      17.079      2.953   5.784  0.00117 ** 
.lin  126.033      7.173  17.570 2.18e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.25 on 6 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 2.143e-06

4.3.2 A pedestrian approach

## Fig. 4.1.
data(segreg)
plot(C ~ Temp, data = segreg, 
     xlab = Mean ~ temperature ~ (degree ~ F),
     ylab = "Energy consumption (kWh)")

profRSS1 <- function(gamma) {
  deviance(lm(C ~ pmax(0, Temp - gamma), data = segreg))
}
profRSS2 <- Vectorize(profRSS1, "gamma")

l <- lm( C ~ pmax(0, Temp - 40), data = segreg)
summary(l)

Call:
lm(formula = C ~ pmax(0, Temp - 40), data = segreg)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.472  -2.155  -0.148   2.678  11.775 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        74.43931    1.16428  63.936  < 2e-16 ***
pmax(0, Temp - 40)  0.53635    0.06226   8.614 2.27e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.312 on 37 degrees of freedom
Multiple R-squared:  0.6673,    Adjusted R-squared:  0.6583 
F-statistic: 74.21 on 1 and 37 DF,  p-value: 2.268e-10

plot(segreg$Temp, segreg$C)
lines(segreg$Temp, predict(l), col = 2)

## Fig. 4.2. 
plot(profRSS2(Temp) ~ Temp, 
     data = segreg,
     type = "l", 
     xlab = expression(gamma),
     ylab = "Profile RSS")

## 좀더 부드럽게 만듬
plot(profRSS2(Temp) ~ Temp, 
     data = data.frame(Temp = seq(min(segreg$Temp), 
                                  max(segreg$Temp), 
                                  length.out = 1001)),
     type = "l", 
     xlab = expression(gamma),
     ylab = "Profile RSS")

4.4 Fitting models with several predictor variables

4.4.1 Two-dimensional predictor

RScompetition.m1 <- nls(
  biomass ~ a/(1 + b * (x + c * z)),
  data = RScompetition,
  start = list(a = 20, b = 1, c = 1)
)
summary(RScompetition.m1)

Formula: biomass ~ a/(1 + b * (x + c * z))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  25.9144     2.3393  11.078 1.42e-14 ***
b   0.1776     0.0355   5.004 8.68e-06 ***
c   1.1349     0.2964   3.829 0.000387 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.673 on 46 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.437e-06

virDensity <- with(RScompetition, x + coef(RScompetition.m1)[3] * z)
virDenVal <- seq(0, max(virDensity), length.out = 100)
biomassVal <- predict(RScompetition.m1, data.frame(x = virDenVal, z = 0))

## Fig. 4.3. 
plot(biomassVal ~ virDenVal, type = "l", 
     ylab = "Biomass of sensitive biotype (g/plant)",
     xlab = Virtual ~ density ~ (plants/m^2))
with(RScompetition, points(biomass ~ virDensity))

4.4.2 General least-squares minimisation

## Fig. 4.4. 
par(pty = "s")
plot(Q ~ I, data = IQsig)
theta <- 0:360 * (pi/180)
lines(cos(theta), sin(theta))

IQsig.m1 <- nls(
  ~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) + 
       gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2, 
  data = IQsig, 
  start = list(I0 = -0.005, gamma = 1, phi = -0.005, Q0 = -0.005, rho = 1))
summary(IQsig.m1)

Formula: 0 ~ ((I - I0)^2 - 2 * gamma * sin(phi) * (I - I0) * (Q - Q0) + 
    gamma * gamma * (Q - Q0)^2) - (rho * gamma * cos(phi))^2

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
I0    -0.002259   0.002058  -1.097    0.278    
gamma  1.007053   0.003719 270.764  < 2e-16 ***
phi   -0.053540   0.004977 -10.758 5.02e-14 ***
Q0    -0.002481   0.002278  -1.089    0.282    
rho    0.974295   0.002540 383.553  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01917 on 45 degrees of freedom

Number of iterations to convergence: 3 
Achieved convergence tolerance: 4.699e-07

5 Model Diagnostics

5.2 Checking the mean structure

5.2.1 Plot of the fitted regression curve

## Fig. 5.1.
plot(p ~ T, data = vapCO, log = "y",
     xlab = "Temperature (K)",
     ylab = "Pressure (Pa)")

vapCO.m1 <- nls(
  log(p) ~ A - B/(C + T), 
  data = vapCO,
  start = list(A = 10, B = 100, C = -10))
summary(vapCO.m1)

Formula: log(p) ~ A - B/(C + T)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
A  19.0923     0.2312   82.56  < 2e-16 ***
B 497.8246    27.7002   17.97 4.84e-10 ***
C -16.1516     1.5437  -10.46 2.19e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09951 on 12 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 3.585e-06

## Fig. 5.2.
plot(p ~ T, data = vapCO, log = "y",
     xlab = "Temperature (K)",
     ylab = "Pressure (Pa)")
lines(vapCO$T, exp(fitted(vapCO.m1)), col = 2)

## Fig. 5.3. 
plot(weight ~ conc, data = lettuce,
     xlab = "Concentration (mg/l)",
     ylab = "Biomass (g)", log = "x")

Warning in xy.coords(x, y, xlabel, ylabel, log) :
  2 x values <= 0 omitted from logarithmic plot

5.2.2 Residual plots

## Fig. 5.4.
plot(fitted(vapCO.m1), residuals(vapCO.m1),
     xlab = "Fitted Values", ylab = "Residuals")
abline(a = 0, b = 0)

## Fig. 5.5. 
plot(rootl ~ conc, data = ryegrass,
     xlab = "Concentration (mM)",
     ylab = "Root length (cm)")

ryegrass.m1 <- lm(rootl ~ as.factor(conc), data = ryegrass)
summary(ryegrass.m1)

Call:
lm(formula = rootl ~ as.factor(conc), data = ryegrass)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.08968 -0.17121  0.01577  0.16847  1.21865 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          7.74935    0.22535  34.388  < 2e-16 ***
as.factor(conc)0.94 -0.07606    0.39032  -0.195  0.84780    
as.factor(conc)1.88 -1.33479    0.39032  -3.420  0.00327 ** 
as.factor(conc)3.75 -4.73466    0.39032 -12.130 8.53e-10 ***
as.factor(conc)7.5  -6.71542    0.39032 -17.205 3.45e-12 ***
as.factor(conc)15   -7.07018    0.39032 -18.114 1.50e-12 ***
as.factor(conc)30   -7.44601    0.39032 -19.077 6.47e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.552 on 17 degrees of freedom
Multiple R-squared:  0.9791,    Adjusted R-squared:  0.9718 
F-statistic:   133 on 6 and 17 DF,  p-value: 2.488e-13

ryegrass.m2 <- nls(
  rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e)))),
  data = ryegrass,
  start = list(b = 1, c = 0.6, d = 8, e = 3)
)
summary(ryegrass.m2)

Formula: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
b   2.9822     0.4584   6.505 2.43e-06 ***
c   0.4814     0.2095   2.298   0.0325 *  
d   7.7930     0.1897  41.073  < 2e-16 ***
e   3.0580     0.1858  16.456 4.31e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5196 on 20 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 6.33e-07

anova(ryegrass.m2, ryegrass.m1)

Analysis of Variance Table

Model 1: rootl ~ c + (d - c)/(1 + exp(b * +(log(conc) - log(e))))
Model 2: rootl ~ as.factor(conc)
  Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
1     20     5.4002                          
2     17     5.1799  3 0.22035  0.2411 0.8665

Q <- -2 * (logLik(ryegrass.m2) - logLik(ryegrass.m1))
df.Q <- df.residual(ryegrass.m2) - df.residual(ryegrass.m1)
1 - pchisq(Q, df.Q)

'log Lik.' 0.8012878 (df=5)

5.3 Variance homogeneity

5.3.2 Levene’s test

with(ryegrass, leveneTest(rootl, as.factor(conc)))

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  6  1.9266 0.1344
      17               

## Fig. 5.6. 
plot(fitted(vapCO.m1), abs(residuals(vapCO.m1)),
     xlab = "Fitted values", ylab = "Absolute residuals")

5.4 Normal distribution

## Fig. 5.7. 
plot(vapCO.m1)

5.4.1 QQ plot

## Fig. 5.8. 
standardRes <- residuals(ryegrass.m2)/summary(ryegrass.m2)$sigma

par(pty = "s")
qqnorm(standardRes, main = "")
abline(a = 0, b = 1)

5.4.2 Shapiro-Wilk test

shapiro.test(standardRes)

    Shapiro-Wilk normality test

data:  standardRes
W = 0.98232, p-value = 0.9345

5.5 Independence

## Fig. 5.9. 
plot(residuals(vapCO.m1), c(residuals(vapCO.m1)[-1], NA),
     xlab = "Residuals", ylab = "Lagged residuals")

6 Remedies for Model Violations

6.1 Variance modelling

6.1.1 Power-of-the-mean variance model

## Fig. 6.1. 
plot(RGR ~ Day, data = RGRcurve,
     xlab = "Time (days)", ylab = "Relative growth rate (%)")

y0_start = mean(RGRcurve[RGRcurve$Day == 0, "RGR"])
y1 = mean(RGRcurve[RGRcurve$Day == 1, "RGR"])
b_start = 1 / log(y1/y0_start)

RGRcurve.m2 <- gnls(
  RGR ~  y0 * exp(Day/b), 
  data = RGRcurve,
  start = list(y0 = y0_start, b = b_start),
  weights = varPower()
)
summary(RGRcurve.m2)

Generalized nonlinear least squares fit
  Model: RGR ~ y0 * exp(Day/b) 
  Data: RGRcurve 

Variance function:
 Structure: Power of variance covariate
 Formula: ~fitted(.) 
 Parameter estimates:
   power 
1.012821 

Coefficients:

 Correlation: 
  y0   
b 0.797

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.8322566 -0.6356287 -0.1612768  0.6476686  1.8176689 

Residual standard error: 0.1672069 
Degrees of freedom: 41 total; 39 residual

RGRcurve.m2g <- gnls(
  RGR ~  y0 * exp(Day/b), 
  data = RGRcurve,
  start = list(y0 = y0_start, b = b_start)
  # weights = varPower() # without weights
)
summary(RGRcurve.m2g )

Generalized nonlinear least squares fit
  Model: RGR ~ y0 * exp(Day/b) 
  Data: RGRcurve 

Coefficients:

 Correlation: 
  y0   
b 0.962

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-2.4946822 -0.5595546 -0.1065972  0.3415913  2.3870757 

Residual standard error: 0.1009825 
Degrees of freedom: 41 total; 39 residual

AIC(RGRcurve.m2, RGRcurve.m2g )

6.2 Transformations

6.2.1 Transform-both-sides approach

## Fig. 6.2. 
plot(recruits ~ spawners, 
     data = sockeye[-12,],
     xlab = "Number of spawners (thousands)",
     ylab = "Number of recruits (thousands)")

일단 nls를 이용하여 fitting을 한 후

sockeye.m1 <- nls(
  recruits ~ beta1 * spawners * exp(-beta2 * spawners),
  data = sockeye[-12, ],
  start = list(beta1 = 2, beta2 = 0.001)
)
summary(sockeye.m1)

Formula: recruits ~ beta1 * spawners * exp(-beta2 * spawners)

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
beta1 3.9495589  1.0154151   3.890 0.000658 ***
beta2 0.0008526  0.0003570   2.388 0.024812 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 566.1 on 25 degrees of freedom

Number of iterations to convergence: 4 
Achieved convergence tolerance: 3.796e-07

## Fig. 6.3. 
plot(fitted(sockeye.m1), abs(residuals(sockeye.m1)),
     xlab = "Fitted values", ylab = "Absolute residuals")

6.2.2 Finding an appropriate transformation

boxcox.nls를 이용하여 양쪽을 변화시키며 피팅을 한다.

## Fig. 6.4.
sockeye.m2 <- boxcox.nls(sockeye.m1)

bcSummary(sockeye.m2)

Estimated lambda: -0.2 
Confidence interval for lambda: [-0.89, 0.52] 

coef(summary(sockeye.m1))

          Estimate   Std. Error  t value    Pr(>|t|)
beta1 3.9495588725 1.0154150955 3.889600 0.000657508
beta2 0.0008525523 0.0003570088 2.388043 0.024811780

coef(summary(sockeye.m2))

          Estimate   Std. Error  t value     Pr(>|t|)
beta1 3.7767882170 0.6941022768 5.441256 1.195243e-05
beta2 0.0009539343 0.0003061565 3.115839 4.563494e-03

6.3 Sandwich estimators

이 내용은 정확히 이해하지는 못했음.

vcov(sockeye.m1)

             beta1        beta2
beta1 1.0310678161 3.436055e-04
beta2 0.0003436055 1.274553e-07

sandwich(sockeye.m1)

             beta1        beta2
beta1 0.6082739015 2.458975e-04
beta2 0.0002458975 1.162302e-07

coeftest(sockeye.m1)

t test of coefficients:

        Estimate Std. Error t value  Pr(>|t|)    
beta1 3.94955887 1.01541510  3.8896 0.0006575 ***
beta2 0.00085255 0.00035701  2.3880 0.0248118 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

coeftest(sockeye.m1, vcov = sandwich)

t test of coefficients:

        Estimate Std. Error t value  Pr(>|t|)    
beta1 3.94955887 0.77991916  5.0641 3.158e-05 ***
beta2 0.00085255 0.00034093  2.5007   0.01931 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

6.4 Weighting

6.4.1 Decline in nitrogen content in soil

exp1

## Fig. 6.5. 
plot(Nremaining ~ time, data = exp1,
     xlab = "Time (years)", ylab = "Nitroge content (%)")

exp1.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1)
exp1.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp1,
               weights = norep/(stdev * stdev))

weights(exp1.m2)

[1] 0.4403928 0.4403928 0.1096776 0.1430179 1.1291355 0.0867694 1.8896447 0.3667661

coef(summary(exp1.m1))

     Estimate Std. Error    t value    Pr(>|t|)
a1 49.6945735  6.1714388  8.0523481 0.001291449
a2  0.2182533  0.3162611  0.6901049 0.528080839
b1 32.2059455  6.1631160  5.2255946 0.006402929
b2 -3.1242900  0.4881789 -6.3998879 0.003061147

coef(summary(exp1.m2))

      Estimate Std. Error     t value    Pr(>|t|)
a1 48.95580529  5.8454395  8.37504275 0.001111764
a2 -0.02288415  0.2953566 -0.07747973 0.941962763
b1 30.67388498  6.4550371  4.75193009 0.008958175
b2 -3.15057479  0.5205589 -6.05229290 0.003760909

exp2

exp2.m1 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2)
exp2.m2 <- nls(Nremaining ~ SSbiexp(time, a1, a2, b1, b2), data = exp2,
               weights = norep/(stdev^2))

## Fig. 6.6. 
plot(Nremaining ~ time, data = exp2,
     xlab = "Time (years)", ylab = "Nitrogen content (%)")
timeVal <- with(exp2, seq(min(time), max(time), length.out = 101))
lines(timeVal, predict(exp2.m1, newdata = data.frame(time = timeVal)), 
      lty = 2, col = 2)

lines(timeVal, predict(exp2.m2, newdata = data.frame(time = timeVal)), 
      lty = 3, col = 3)

coef(summary(exp2.m1))

     Estimate Std. Error     t value    Pr(>|t|)
a1 37.1828528  7.3946625   5.0283367 0.007342087
a2  0.3984263  0.4676840   0.8519136 0.442262753
b1 56.9940958  7.2995778   7.8078620 0.001451984
b2 -2.4593036  0.2429193 -10.1239548 0.000535818

coef(summary(exp2.m2))

     Estimate Std. Error     t value     Pr(>|t|)
a1 28.9158392  8.3745109   3.4528392 2.598700e-02
a2  0.3128281  0.8122243   0.3851499 7.197304e-01
b1 59.7777230  3.7595533  15.9002197 9.144765e-05
b2 -2.3567707  0.1311759 -17.9664898 5.641332e-05

summary(exp1.m1)

Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a1  49.6946     6.1714   8.052  0.00129 **
a2   0.2183     0.3163   0.690  0.52808   
b1  32.2059     6.1631   5.226  0.00640 **
b2  -3.1243     0.4882  -6.400  0.00306 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.844 on 4 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 5.364e-06

summary(exp1.m2)

Formula: Nremaining ~ SSbiexp(time, a1, a2, b1, b2)

Parameters:
   Estimate Std. Error t value Pr(>|t|)   
a1 48.95581    5.84544   8.375  0.00111 **
a2 -0.02288    0.29536  -0.077  0.94196   
b1 30.67388    6.45504   4.752  0.00896 **
b2 -3.15057    0.52056  -6.052  0.00376 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.785 on 4 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 9.037e-06

7 Uncertainty, Hypothesis Testing, and Model Selection

7.1 Profile likelihood

L.minor.m1 <- update(L.minor.m1, trace = FALSE)
L.minor.m1pro <- profile(L.minor.m1)

## Fig. 7.1. 
plot(L.minor.m1pro)

## Fig. 7.2. 
plot(L.minor.m1pro, absVal = FALSE)

confint(L.minor.m1)

Waiting for profiling to be done...

        2.5%     97.5%
K   11.45142  24.86684
Vm 110.64232 143.75648

7.2 Bootstrap

set.seed(222)
L.minor.m1boot <- nlsBoot(L.minor.m1)
summary(L.minor.m1boot)

------
Bootstrap statistics
    Estimate Std. error
K   17.37534   2.730157
Vm 126.64954   6.522141

------
Median of bootstrap estimates and percentile confidence intervals
      Median      2.5%     97.5%
K   17.25307  12.30746  23.20204
Vm 126.08120 115.00875 139.46932

## Fig. 7.3. 
par(mfrow = c(1, 2))
qqnorm(L.minor.m1boot$coefboot[, 1], main = "K")
qqnorm(L.minor.m1boot$coefboot[, 2], main = "Vm")

7.3 Wald confidence intervals

confint2(L.minor.m1)

        2.5 %    97.5 %
K    9.853828  24.30416
Vm 108.480813 143.58470

confint2(L.minor.m1, level = 0.99)

       0.5 %    99.5 %
K   6.131815  28.02617
Vm 99.439005 152.62651

7.4 Estimating derived parameters

deltaMethod(L.minor.m1, "Vm/(4*K)")

           Estimate     SE  2.5 % 97.5 %
Vm/(4 * K)   1.8449 0.2363 1.3817  2.308

7.5 Nested models

secalonic

## Fig. 7.4.
plot(rootl ~ dose, data = secalonic,
     xlab = "Dose (mM)", ylab = "Root length (cm)")

7.5.1 Using t-tests

secalonic.m1 <- nls(rootl ~ SSfpl(dose, a, b, c, d), data = secalonic)
summary(secalonic.m1)

Formula: rootl ~ SSfpl(dose, a, b, c, d)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a 6.053612   0.395467  15.308 0.000606 ***
b 0.353944   0.194089   1.824 0.165722    
c 0.075188   0.005911  12.721 0.001048 ** 
d 0.029350   0.006621   4.433 0.021333 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2028 on 3 degrees of freedom

Number of iterations to convergence: 0 
Achieved convergence tolerance: 8.102e-06

7.5.2 Using F-tests

secalonic.m2 <- nls(rootl ~ SSlogis(dose, a, c, d), data = secalonic)
summary(secalonic.m2)

Formula: rootl ~ SSlogis(dose, a, c, d)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
a  6.344853   0.564690  11.236 0.000357 ***
c  0.076959   0.009397   8.190 0.001211 ** 
d -0.035740   0.007244  -4.933 0.007854 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2546 on 4 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.999e-07

anova(secalonic.m2, secalonic.m1)

Analysis of Variance Table

Model 1: rootl ~ SSlogis(dose, a, c, d)
Model 2: rootl ~ SSfpl(dose, a, b, c, d)
  Res.Df Res.Sum Sq Df  Sum Sq F value Pr(>F)
1      4    0.25924                          
2      3    0.12341  1 0.13582  3.3016 0.1668

7.6 Non-nested models

M.merluccius

plot(num.fish ~ spawn.biomass, data = M.merluccius)

4가지 조금씩 다른 모델을 비교

M.merluccius.bh <- nls(
  num.fish ~ spawn.biomass * alpha/(1 + spawn.biomass/k),
  data = M.merluccius, start = list(alpha = 5, k = 50)
)

M.merluccius.de <- nls(
  num.fish ~ spawn.biomass * alpha * (1 - c * spawn.biomass/k)^(1/c),
  data = M.merluccius, start = list(alpha = 4.4, k = 106, c = 0.86)
)

M.merluccius.ri <- nls(
  num.fish ~ spawn.biomass * alpha * exp(-spawn.biomass/k),
  data = M.merluccius, start = list(alpha = 5, k = 50)
)

M.merluccius.sh <- nls(
  num.fish ~ spawn.biomass * alpha/(1 + (spawn.biomass/k)^c),
  data = M.merluccius, start = list(alpha = 3.87, k = 61.72, c = 2.25),
  control = nls.control(maxiter = 100)
)

Residual standard error를 계산하면,

summary(M.merluccius.bh)$sigma

[1] 14.69956

summary(M.merluccius.de)$sigma

[1] 14.78346

summary(M.merluccius.ri)$sigma

[1] 14.30425

summary(M.merluccius.sh)$sigma

[1] 14.56758

AIC를 비교해보면,

AIC(M.merluccius.bh)

[1] 127.0562

AIC(M.merluccius.de)

[1] 128.0263

AIC(M.merluccius.ri)

[1] 126.2383

AIC(M.merluccius.sh)

[1] 127.585

Burnham and Anderson (2002, pp. 70–72)은 한 모델이 다른 모델에 비하여 절대적으로 우위에 있기 위해서는 AIC의 차이가 10 이상의 차이가 필요하다는 경험적 법칙을 제시하였다. 따라서 4 모델 사이에서 큰 차이는 없다.

8 Grouped Data

8.1 Fitting grouped data models

Puromycin

## Fig. 8.1.
xyplot(rate ~ conc | state, data = Puromycin,
       xlab = "Substrate concentration (ppm)",
       ylab = "Reaction rates\n(counts/min/min)")

8.1.1 Using nls()

Puromycin.m1 <- nls(
  rate ~ Vm[state] * conc/(K[state] + conc),
  data = Puromycin,
  start = list(K = c(0.1, 0.1), Vm = c(200, 200))
)
summary(Puromycin.m1)

Formula: rate ~ Vm[state] * conc/(K[state] + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K1  6.412e-02  7.877e-03   8.141 1.29e-07 ***
K2  4.771e-02  8.281e-03   5.761 1.50e-05 ***
Vm1 2.127e+02  6.608e+00  32.185  < 2e-16 ***
Vm2 1.603e+02  6.896e+00  23.242 2.04e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.4 on 19 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.568e-06

8.1.2 Using gnls()

Puromycin.m2 <- gnls(
  rate ~ Vm * conc/(K + conc),
  data = Puromycin,
  start = list(Vm = c(200, 200), K = c(0.1, 0.1)),
  params = list(Vm ~ state - 1, K ~ state - 1)
)
summary(Puromycin.m2)

Generalized nonlinear least squares fit
  Model: rate ~ Vm * conc/(K + conc) 
  Data: Puromycin 

Coefficients:

 Correlation: 
                  Vm.sttt Vm.sttn K.sttt
Vm.stateuntreated 0.000                 
K.statetreated    0.765   0.000         
K.stateuntreated  0.000   0.777   0.000 

Standardized residuals:
        Min          Q1         Med          Q3         Max 
-1.32632770 -0.52742178 -0.06890298  0.51295804  2.44557119 

Residual standard error: 10.40003 
Degrees of freedom: 23 total; 19 residual

8.1.3 Using nlsList()

Puromycin.m3 <- nlsList(
  rate ~ SSmicmen(conc, Vm, K) | state,
  data = Puromycin
)
summary(Puromycin.m3)

Call:
  Model: rate ~ SSmicmen(conc, Vm, K) | state 
   Data: Puromycin 

Coefficients:
   Vm 
          Estimate Std. Error  t value     Pr(>|t|)
treated   212.6837   6.608092 32.18534 3.241160e-11
untreated 160.2801   6.896019 23.24241 1.384618e-09
   K 
            Estimate  Std. Error  t value     Pr(>|t|)
treated   0.06412123 0.007876786 8.140532 1.565136e-05
untreated 0.04770829 0.008281171 5.761056 1.727047e-04

Residual standard error: 10.40003 on 19 degrees of freedom

Puromycin2 <- groupedData(rate ~ conc | state, data = Puromycin)
Puromycin.m4 <- nlsList(
  rate ~ SSmicmen(conc, a, b), data = Puromycin2
)
summary(Puromycin.m4)

Call:
  Model: rate ~ SSmicmen(conc, a, b) | state 
   Data: Puromycin2 

Coefficients:
   a 
          Estimate Std. Error  t value     Pr(>|t|)
untreated 160.2801   6.896019 23.24241 1.384618e-09
treated   212.6837   6.608092 32.18534 3.241160e-11
   b 
            Estimate  Std. Error  t value     Pr(>|t|)
untreated 0.04770829 0.008281171 5.761056 1.727047e-04
treated   0.06412123 0.007876786 8.140532 1.565136e-05

Residual standard error: 10.40003 on 19 degrees of freedom

8.2 Model reduction and parameter models

8.2.1 Comparison of entire groups

Puromycin.m5 <- nls(
  rate ~ Vm * conc/(K + conc), data = Puromycin,
  start = list(K = 0.1, Vm = 200))
summary(Puromycin.m5)

Formula: rate ~ Vm * conc/(K + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K    0.06039    0.01077   5.608 1.45e-05 ***
Vm 190.80620    8.76458  21.770 6.84e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 18.61 on 21 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 6.259e-06

anova(Puromycin.m5, Puromycin.m1)

Analysis of Variance Table

Model 1: rate ~ Vm * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value    Pr(>F)    
1     21     7276.5                                
2     19     2055.1  2 5221.5  24.138 6.075e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

AIC(Puromycin.m5, Puromycin.m1)

8.2.2 Comparison of specific parameters

Puromycin.m6 <- nls(
  rate ~ Vm * conc/(K[state] + conc), 
  data = Puromycin, 
  start = list(K = c(0.1, 0.1), Vm = 200))
summary(Puromycin.m6)

Formula: rate ~ Vm * conc/(K[state] + conc)

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K1   0.04801    0.00874   5.493 2.24e-05 ***
K2   0.08865    0.01448   6.122 5.54e-06 ***
Vm 194.10328    7.40464  26.214  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.52 on 20 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 8.957e-06

anova(Puromycin.m6, Puromycin.m1)

Analysis of Variance Table

Model 1: rate ~ Vm * conc/(K[state] + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value    Pr(>F)    
1     20     4815.4                                
2     19     2055.1  1 2760.3  25.521 7.082e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Puromycin.m7 <- nls(
  rate ~ Vm[state] * conc/(K + conc), 
  data = Puromycin, start = list(K = 0.1, Vm = c(200, 200))
)
summary(Puromycin.m7)

Formula: rate ~ Vm[state] * conc/(K + conc)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
K     0.05797    0.00591   9.809 4.37e-09 ***
Vm1 208.63004    5.80399  35.946  < 2e-16 ***
Vm2 166.60408    5.80743  28.688  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.59 on 20 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 1.364e-06

anova(Puromycin.m7, Puromycin.m1)

Analysis of Variance Table

Model 1: rate ~ Vm[state] * conc/(K + conc)
Model 2: rate ~ Vm[state] * conc/(K[state] + conc)
  Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1     20     2240.9                         
2     19     2055.1  1 185.84  1.7182 0.2056

AIC로 비교해 보아도 Puromycin.m7이 더욱 좋다.

AIC(Puromycin.m1, Puromycin.m5, Puromycin.m6, Puromycin.m7)

8.3 Common control

G.aparine

## Fig. 8.2. 
xyplot(drymatter ~ dose | as.factor(treatment),
       data = G.aparine,
       xlab = "Dose (g/ha)", ylab = "Dry matter (mg/pot)")

이미 데이터에서 treatment 0가 없어졌음. 아마도 데이터를 업데이트 한 것으로 추측됨.

G.aparine$treatment2 <- factor(G.aparine$treatment)
#levels(G.aparine$treatment2) <- c("0", "1", "2")
G.aparine

G.aparine.m1 <- nls(
  drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] * (log(dose) - log(e[treatment2])))),
  data = G.aparine, 
  start = list(b = c(2, 2), c = c(500, 100), d = 1000, e = c(50, 100)))
summary(G.aparine.m1)

Formula: drymatter ~ c[treatment2] + (d - c[treatment2])/(1 + exp(b[treatment2] * 
    (log(dose) - log(e[treatment2]))))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b1   1.6129     0.3462   4.659 5.33e-06 ***
b2   1.7510     0.2303   7.603 7.07e-13 ***
c1 509.5035    23.4379  21.738  < 2e-16 ***
c2 151.9164    27.2114   5.583 6.56e-08 ***
d  984.8883    12.7957  76.970  < 2e-16 ***
e1  50.8000     7.8465   6.474 5.58e-10 ***
e2  93.4466     8.1626  11.448  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 111.4 on 233 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 9.945e-06

8.4 Prediction

concValues <- with(Puromycin, seq(min(conc), max(conc, length.out = 11)))
concValues

 [1]  0.02  1.02  2.02  3.02  4.02  5.02  6.02  7.02  8.02  9.02 10.02

stateVal1 <- levels(Puromycin$state)
stateVal1

[1] "treated"   "untreated"

csValues1 <- expand.grid(conc = concValues, state = stateVal1)
csValues1

predict(Puromycin.m7, newdata = csValues1)

 [1]  53.51423 197.41021 202.80963 204.70062 205.66419 206.24825 206.64012
 [8] 206.92127 207.13280 207.29773 207.42993  42.73445 157.64435 161.95611
[15] 163.46619 164.23566 164.70207 165.01500 165.23951 165.40844 165.54015
[22] 165.64572

stateVal2 <- factor("untreated", levels = c("treated", "untreated"))
stateVal2

[1] untreated
Levels: treated untreated

csValues2 <- data.frame(conc = concValues, state = stateVal2)
csValues2

predict(Puromycin.m1, newdata = csValues2)

 [1]  47.3444 153.1183 156.5819 157.7874 158.4002 158.7711 159.0198 159.1981
 [9] 159.3322 159.4367 159.5205

8.5 Nonlinear mixed models

vinclozolin

## Fig. 8.3.
xyplot(effect ~ conc | exper, data = vinclozolin,
       xlab = expression(paste("Concentration (", mu, "M)")),
       ylab = "Luminescence (LU")

LL3.formula <- effect ~ d/(1 + exp(b * (log(conc) - log(e))))

vinclozolin.e1.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10509,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e2.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10821,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e3.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10828,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e4.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 10904,  
  start = list(b = 1, d = 2700, e = 0.03))
vinclozolin.e5.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 11023,  
  start = list(b = 1, d = 1000, e = 0.26))
vinclozolin.e6.m <- nls(
  LL3.formula, data = vinclozolin, subset = exper == 11106,  
  start = list(b = 0.5, d = 2600, e = 0.02))

## Fig. 8.4. 
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
     pch = as.numeric(exper), 
     log = "x",
     xlim = c(1e-04, 10), 
     xlab = expression(paste("Concentration(", mu, "M)")), 
     ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)

lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)

lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)

lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")

vinclozolin.m1 <- nlme(
  effect ~ d/(1 + exp(b * (log(conc) - log(e)))), 
  data = vinclozolin,
  fixed = list(b ~ 1, d ~ 1, e ~ 1), 
  random = d ~ 1 | exper, 
  start = c(1, 1000, 0.1))
summary(vinclozolin.m1)

Nonlinear mixed-effects model fit by maximum likelihood
  Model: effect ~ d/(1 + exp(b * (log(conc) - log(e)))) 
  Data: vinclozolin 

Random effects:
 Formula: d ~ 1 | exper
               d Residual
StdDev: 695.3775 150.9623

Fixed effects:  list(b ~ 1, d ~ 1, e ~ 1) 

 Correlation: 
  b      d     
d -0.051       
e  0.529 -0.145

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.83842036 -0.58871949 -0.06161384  0.54652838  2.07902663 

Number of Observations: 53
Number of Groups: 6 

## Fig. 8.5.
#par(pty = "s")
plot(effect ~ conc, data = vinclozolin,
     pch = as.numeric(exper), 
     log = "x",
     xlim = c(1e-04, 10), 
     xlab = expression(paste("Concentration(", mu, "M)")), 
     ylab = "Luminescence(LU)")
concVec <- exp(seq(log(1e-04), log(10), length.out = 50))
lines(concVec, predict(vinclozolin.e1.m, data.frame(conc = concVec)), lty = 2)

lines(concVec, predict(vinclozolin.e2.m, data.frame(conc = concVec)), lty = 3)
lines(concVec, predict(vinclozolin.e3.m, data.frame(conc = concVec)), lty = 4)

lines(concVec, predict(vinclozolin.e4.m, data.frame(conc = concVec)), lty = 5)
lines(concVec, predict(vinclozolin.e5.m, data.frame(conc = concVec)), lty = 6)

lines(concVec, predict(vinclozolin.e6.m, data.frame(conc = concVec)), lty = "3313")
lines(concVec, 
      predict(vinclozolin.m1, newdata = data.frame(conc = concVec), level = 0), 
      lty = 1, lwd = 3)

—-